Why the Z-Transform?

• The Fourier Transform assumes signals are stable and well-behaved
• But some biosignals or systems may not be absolutely summable
• The Z-Transform generalizes the Fourier Transform
• Useful for analyzing discrete-time systems, especially when stability and causality matter

Definition

Let \(x[n]\) be a discrete-time signal.

The Z-Transform is defined as:

\[X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}\]

Where: - \(z \in \mathbb{C}\) is a complex variable - \(z = re^{j\omega}\)

Region of Convergence (ROC)

• The Z-Transform converges only for certain values of \(z\)
• The set of \(z\) for which the series converges is the ROC
• ROC is critical for system stability and causality

Z-Plane Representation

• Poles: values of \(z\) where \(X(z) \to \infty\)
• Zeros: values where \(X(z) = 0\)
• Visualization of poles and zeros helps in understanding system behavior

(-1.5, 1.5)

(-1.5, 1.5)

Relationship with Fourier Transform

If the ROC includes the unit circle, \(|z| = 1\), then:

\[X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}\]

So the Fourier Transform is a special case of the Z-Transform.

Example Signal

Let:

\[x[n] = a^n u[n]\]

Where: - \(a \in \mathbb{R}\) - \(u[n]\) is the unit step function (0 for \(n<0\), 1 for \(n\geq 0\))

Step 1: Z-Transform

Apply the definition:

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \sum_{n=0}^{\infty} (az^{-1})^n\]

This is a geometric series:

\[X(z) = \frac{1}{1 - az^{-1}} = \frac{z}{z - a}\]

What is a Geometric Series?

A geometric series is a sum of terms where each term is multiplied by the same constant:

\[S = \sum_{n=0}^{\infty} r^n = 1 + r + r^2 + r^3 + \cdots\]

The value of \(r\) determines whether this sum converges (has a finite limit) or diverges.

Goal

Understand why this series:

\[\sum_{n=0}^{\infty} r^n\]

converges if and only if:

\[\boxed{|r| < 1}\]

Partial Sums

Let’s consider the sum up to the \(N\)-th term:

\[S_N = \sum_{n=0}^{N} r^n = 1 + r + r^2 + \cdots + r^N\]

This has a closed-form expression:

\[S_N = \frac{1 - r^{N+1}}{1 - r} \quad \text{for } r \neq 1\]

Taking the Limit

To find the sum of the infinite series, take the limit as \(N \to \infty\):

\[S = \lim_{N \to \infty} S_N = \lim_{N \to \infty} \frac{1 - r^{N+1}}{1 - r}\]

Case 1: \(|r| < 1\)

If \(|r| < 1\), then:

\[r^{N+1} \to 0 \quad \text{as } N \to \infty\]

So the sum becomes:

\[S = \frac{1}{1 - r}\]

✅ The geometric series converges.

Case 2: \(|r| \geq 1\)

• If \(r = 1\), then \(S_N = N + 1 \to \infty\)
• If \(r = -1\), the sum oscillates: \(1 - 1 + 1 - 1 + \cdots\)
• If \(|r| > 1\), then \(r^{N+1} \to \infty\)

❌ In all cases: the series diverges

Intuition

• When \(|r| < 1\), each term \(r^n\) gets smaller and smaller
• Their total sum settles to a finite number
• When \(|r| \geq 1\), the terms don’t vanish — the sum keeps growing or oscillating

Example: \(r = 0.5\)

\[S = 1 + 0.5 + 0.25 + 0.125 + \cdots = \frac{1}{1 - 0.5} = 2\]

Every term adds less. The sum “flattens out.”

Why This Matters

The Z-Transform often gives us geometric series like:

\[\sum_{n=0}^{\infty} (az^{-1})^n\]

This converges only if:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

So, understanding convergence of geometric series = understanding ROC in Z-transforms.

Summary

\[\sum_{n=0}^{\infty} r^n = \frac{1}{1 - r} \quad \text{if } |r| < 1\]

Step 2: Region of Convergence (ROC)

For convergence of the geometric series:

\[|az^{-1}| < 1 \Rightarrow |z| > |a|\]

Therefore, the ROC is:

\[\boxed{|z| > |a|}\]

Interpretation

• Causal signal (defined for \(n \geq 0\))
• ROC is outside the outermost pole
• Stable system only if ROC includes the unit circle \(|z| = 1\)

Example: \(a = 0.5\)

\[x[n] = (0.5)^n u[n]\]

Z-Transform:

\[X(z) = \sum_{n=0}^{\infty} (0.5)^n z^{-n} = \frac{z}{z - 0.5}\]

Region of Convergence:

\[\boxed{|z| > 0.5}\]

Summary

Properties of the Z-Transform

• Linearity: \(a x[n] + b y[n] \to aX(z) + bY(z)\)
• Time shifting: \(x[n - k] \to z^{-k} X(z)\)
• Scaling in the z-domain: \(a^n x[n] \to X(z/a)\)
• Convolution: \(x[n] * h[n] \to X(z)H(z)\)

Example

Let \(x[n] = a^n u[n]\), where \(|a| < 1\)

\[X(z) = \sum_{n=0}^{\infty} a^n z^{-n} = \frac{1}{1 - az^{-1}}, \quad \text{ROC: } |z| > |a|\]

Difference Equations in DSP

A difference equation relates input and output values at different time steps.

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Common in: - Digital filters (FIR, IIR) - Signal models in ECG, EEG analysis - Implementation in real-time biosignal systems

Z-Transform of Time-Shifted Terms

The Z-Transform turns time shifts into powers of \(z^{-1}\):

Step 1: Apply Z-Transform

Given:

\[y[n] - a_1 y[n-1] - a_2 y[n-2] = b_0 x[n] + b_1 x[n-1]\]

Apply \(\mathcal{Z} \{ \cdot \}\):

\[Y(z) - a_1 z^{-1} Y(z) - a_2 z^{-2} Y(z) = b_0 X(z) + b_1 z^{-1} X(z)\]

Step 2: Factor and Solve for \(H(z)\)

Group:

\[Y(z)(1 - a_1 z^{-1} - a_2 z^{-2}) = X(z)(b_0 + b_1 z^{-1})\]

Divide both sides:

\[H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1}}{1 - a_1 z^{-1} - a_2 z^{-2}}\]

Example

Given:

\[y[n] - 0.9 y[n-1] = x[n] - 0.5 x[n-1]\]

Z-Transform:

\[Y(z)(1 - 0.9 z^{-1}) = X(z)(1 - 0.5 z^{-1})\]

Transfer Function:

\[H(z) = \frac{1 - 0.5 z^{-1}}{1 - 0.9 z^{-1}}\]

Poles and Zeros

Let’s analyze \(H(z)\):

• Zeros: Roots of the numerator \(\Rightarrow z = 0.5\)
• Poles: Roots of the denominator \(\Rightarrow z = 0.9\)

Practice

Convert this equation:

\[y[n] = 0.6 y[n-1] + x[n] + x[n-1]\]

Find: - \(H(z)\) - Poles and zeros - Plot them in the Z-plane

Application in Biosignal Processing

• Analysis of digital filters for ECG, EEG, etc.
• Design of stable and causal filtering systems
• Useful in difference equation modeling of biosignals

Summary

• Z-Transform is a powerful tool for analyzing discrete systems
• Provides insight into stability, causality, and system behavior
• A generalization of the Fourier Transform
• Crucial in digital signal processing of biosignals
• Z-Transform converts difference equations into algebraic expressions
• Transfer function \(H(z)\) tells us how the system responds to inputs
• Key for digital filter design in biosignal processing

Next Steps

• Practice Z-Transform computations
• Pole-zero plotting exercises
• Application to real biosignal filtering problems

Filter Design

Filter Design

def design_filter(zeros=None, poles=None, gain=1.0):
    """
    Diseña un filtro digital a partir de ceros y/o polos y una ganancia.

    Parámetros:
    - zeros: lista de ceros (raíces del numerador), o None para no incluir
    - poles: lista de polos (raíces del denominador), o None para no incluir
    - gain: ganancia escalar del filtro

    Devuelve:
    - b: coeficientes del numerador
    - a: coeficientes del denominador
    """
    # Si no se pasan ceros, asumimos un FIR trivial (b = [gain])
    if zeros:
        b = gain * np.poly(zeros)
    else:
        b = np.array([gain], dtype=float)

    # Si no se pasan polos, asumimos sistema FIR (a = [1])
    if poles:
        a = np.poly(poles)
    else:
        a = np.array([1.0], dtype=float)

    return b, a

Filter Design

def gaussian(x, mu, sigma, A):
    """
    Genera una función gaussiana.

    Parámetros:
    - x: array de tiempos
    - mu: posición central de la gaussiana
    - sigma: desviación estándar
    - A: amplitud
    """
    return A * np.exp(-((x - mu) ** 2) / (2 * sigma**2))

Filter Design

def simulate_ecg(duration=10.0, fs=500, heart_rate=60):
    """
    Simula un ECG sintético basado en la superposición de ondas gaussianas.

    Parámetros:
    - duration: duración de la señal en segundos
    - fs: frecuencia de muestreo en Hz
    - heart_rate: frecuencia cardiaca en latidos por minuto

    Devuelve:
    - t: vector de tiempos
    - ecg: señal simulada de ECG en mV
    """
    dt = 1 / fs
    t = np.arange(0, duration, dt)
    rr = 60 / heart_rate  # intervalo RR en segundos

    # Inicializar señal
    ecg = np.zeros_like(t)

    # Parámetros de las ondas (posiciones relativas en segundos)
    # P wave
    p_amp, p_dur, p_delay = 0.25, 0.09, 0.16
    # Q wave
    q_amp, q_dur, q_delay = -0.05, 0.066, 0.166
    # R wave
    r_amp, r_dur, r_delay = 1.6, 0.1, 0.166
    # S wave
    s_amp, s_dur, s_delay = -0.25, 0.066, 0.19
    # T wave
    t_amp, t_dur, t_delay = 0.35, 0.142, 0.36

    # Generar cada latido
    for beat_start in np.arange(0, duration, rr):
        mask = (t >= beat_start) & (t < beat_start + rr)
        tb = t[mask] - beat_start
        ecg[mask] += gaussian(tb, p_delay, p_dur / 2, p_amp)
        ecg[mask] += gaussian(tb, q_delay, q_dur / 2, q_amp)
        ecg[mask] += gaussian(tb, r_delay, r_dur / 2, r_amp)
        ecg[mask] += gaussian(tb, s_delay, s_dur / 2, s_amp)
        ecg[mask] += gaussian(tb, t_delay, t_dur / 2, t_amp)

    return t, ecg+0.1*np.cos(2*np.pi*60*t)

Filter Design

# Parámetros de simulación
DURATION = 10    # segundos
FS = 500         # Hz
HR = 70          # latidos por minuto

# Generar señal
t, ecg_signal = simulate_ecg(duration=DURATION, fs=FS, heart_rate=HR)

# Graficar resultado
plt.figure(figsize=(12, 4))
plt.plot(t, ecg_signal, linewidth=1)
plt.title(f'Señal de ECG sintética ({HR} bpm)')
plt.xlabel('Tiempo (s)')
plt.ylabel('Amplitud (mV)')
plt.grid(True)
plt.tight_layout()
plt.show()

Filter Design

n = len(ecg_signal)
yf = np.fft.fft(ecg_signal)
xf = np.fft.fftfreq(n, 1/fs)[:n//2]
plt.plot(xf, 2.0/n * np.abs(yf[0:n//2]))
plt.grid()
plt.show()

Filter Design

fc = 60
b,a =design_filter(zeros=[np.exp(1j*2*np.pi*fc/fs), np.exp(-1j*2*np.pi*fc/fs)])
print(a)
print(b)

w, h = sig.freqz(a, b, worN=8000)
plt.plot(w/np.pi*fs/2, 20*np.log10(abs(h)))
plt.grid()
plt.xlabel('Frequency (Hz)')

Filter Design

Filter Design

You obtain \(h_d[n]\) as the inverse discrete–time Fourier transform of your ideal frequency response \(H_d(e^{j\omega})\). For a low-pass filter with cutoff \(\omega_c\),

\[ H_d(e^{j\omega}) = \begin{cases} 1, & |\omega|\le\omega_c,\\ 0, & \text{otherwise}. \end{cases} \]

By definition of the inverse DTFT,

\[ h_d[n] \;=\; \frac{1}{2\pi}\int_{-\pi}^{\pi} H_d(e^{j\omega})\,e^{j\omega n}\,d\omega \;=\;\frac{1}{2\pi}\int_{-\omega_c}^{\omega_c} e^{j\omega n}\,d\omega. \]

Relating to sampling frequency

If your cutoff is specified in Hz, \(f_c\), and sampling rate is \(f_s\), then

\[ \omega_c = 2\pi\,\frac{f_c}{f_s}, \]

so you can write

\[ h_d[n] =\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{\pi\,n} =\;2\;\frac{f_c}{f_s}\;\frac{\sin\!\bigl(2\pi \frac{f_c}{f_s}\,n\bigr)}{2\pi \frac{f_c}{f_s}\,n} =2\frac{f_c}{f_s}\,\mathrm{sinc}\!\Bigl(2\frac{f_c}{f_s}\,n\Bigr), \]

where we define the normalized sinc as \(\mathrm{sinc}(x)=\frac{\sin(\pi x)}{\pi x}\).

Key takeaway

• \(h_d[n]\) is exactly the inverse‐DTFT of the ideal (“brick‐wall”) frequency specification.
• It produces a sinc-shaped impulse response of infinite length.
• Truncation (with a window) makes it realizable as a finite-length FIR.

Filter Design

The typical characteristics of classic windows are:

• Rectangular: narrowest main lobe (width ≈ 4π/M rad) but highest sidelobes (first sidelobe ≈ –13 dB, stop-band attenuation ≈ 21 dB). It gives the steepest transition for a given filter order, at the expense of poorer stop-band rejection.

• Bartlett (triangular): somewhat wider main lobe (≈ 8π/M), sidelobes ≈ –25 dB.

• Hann: main lobe ≈ 8π/M, sidelobes ≈ –31 dB (better rejection than rectangular but smoother transitions).

• Hamming: main lobe ≈ 8π/M, sidelobes ≈ –41 dB (minimum stop-band attenuation ≈ 53 dB). Very popular for its good compromise between transition width and stop-band attenuation.

• Blackman: wider main lobe (≈ 12π/M) but very low sidelobes (≈ –57 dB, attenuation ≈ 74 dB).

• Kaiser: allows selection of a parameter β to control sidelobe attenuation, offering flexibility. Approximately, to achieve A dB of attenuation,

  \[ \beta \approx 0.1102\,(A - 8.7)\quad(\text{for }A>50), \]

  and the normalized transition width Δω relates to the order M and β by

  \[ M \approx \frac{A - 8}{2.285\,\Delta\omega} \]

These formulas stem from Kaiser’s approximations and help in sizing the filter.

Windows Forms

(-100.0, 5.0)

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import (
    boxcar,       # Rectangular
    bartlett,    # Triangular
    hann,        # Hann
    hamming,     # Hamming
    blackman,    # Blackman
    kaiser       # Kaiser
)

# Parameters
M = 64               # window length
beta = 8.6           # Kaiser parameter for moderate sidelobe attenuation
nfft = 512           # FFT length for frequency response
fs = 1.0             # normalized sampling rate

# Generate windows
windows = {
    'Rectangular': boxcar(M),
    'Bartlett':    bartlett(M),
    'Hann':        hann(M),
    'Hamming':     hamming(M),
    'Blackman':    blackman(M),
    f'Kaiser (β={beta})': kaiser(M, beta)
}

# Time-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    plt.plot(np.arange(M), w, label=name)
plt.title('Window Functions — Time Domain')
plt.xlabel('Sample index n')
plt.ylabel('Amplitude')
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

# Frequency-domain plot
plt.figure(figsize=(8, 4))
for name, w in windows.items():
    # Compute normalized frequency response
    W = np.fft.fft(w, nfft)
    W_mag = 20 * np.log10(np.abs(W) / np.max(np.abs(W)))
    freqs = np.fft.fftfreq(nfft, d=1/fs)
    # Only plot 0 ≤ f ≤ 0.5 (normalized Nyquist)
    idx = np.logical_and(freqs >= 0, freqs <= 0.5)
    plt.plot(freqs[idx], W_mag[idx], label=name)
plt.title('Window Functions — Frequency Response')
plt.xlabel('Normalized Frequency (cycles/sample)')
plt.ylabel('Magnitude (dB)')
plt.ylim(-100, 5)
plt.legend(loc='upper right')
plt.grid(True)
plt.tight_layout()

plt.show()

Filter Design

N = 101
fc = 30

# 2. Compute the ideal impulse response h_d[n] of a low-pass filter
n = np.arange(N)
M = (N - 1) / 2
# Using the normalized sinc: sinc(x) = sin(pi*x)/(pi*x)
h_d = (2 * fc / FS) * np.sinc(2 * fc * (n - M) / FS)

# 3. Choose a window w[n] (here: Hamming) and truncate
w = np.hamming(N)
h = h_d * w

# 4. Normalize to ensure unity gain at DC
h /= np.sum(h)

# 6. Filter it (only allowed library call)
y = sig.lfilter(h, 1.0, ecg_signal)

# 7. Plot input vs. output
plt.figure(figsize=(8, 4))
plt.plot(t, ecg_signal, label='Original signal')
plt.plot(t, y, label='Filtered signal', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.title('Low-pass FIR (window method) – Cutoff 100 Hz')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

Step 0: The Analog Domain - Laplace Transform

• Context: Continuous-time signals and Linear Time-Invariant (LTI) systems.
• Tool: The Laplace Transform converts differential equations (time domain) to algebraic equations (s-domain).
• Transfer Function \(H(s)\): Represents the system in the s-domain. \[H(s) = \frac{Y(s)}{X(s)}\] where \(s = \sigma + j\Omega\) is the complex frequency (\(\Omega\) = analog angular frequency).

Step 0: The Analog Domain - s-Plane Analysis

• s-Plane: Complex plane where \(H(s)\) is analyzed.
• Poles & Zeros: Roots of the denominator & numerator of \(H(s)\), respectively. Their locations determine filter behavior.
• Frequency Response: Determined by \(H(j\Omega)\) (evaluating \(H(s)\) on the imaginary axis).
• Stability: For a causal system to be stable, all poles must be in the Left-Half Plane (LHP), i.e., \(\text{Re}\{s\} < 0\).

Step 0: The Analog Domain - Filter Prototypes

• Standardized, well-understood analog filter approximations:
  • Butterworth: Maximally flat passband.
  • Chebyshev Type I: Equiripple passband, monotonic stopband.
  • Chebyshev Type II: Monotonic passband, equiripple stopband.
  • Elliptic (Cauer): Equiripple in both passband and stopband (sharpest transition for a given order).
• Design usually starts with a normalized low-pass prototype.

The Bridge: Analog-to-Digital Mapping

• Core Task: Convert the analog filter \(H(s)\) (s-plane) into a digital filter \(H(z)\) (z-plane).
• Need a mapping that relates the continuous frequency \(\Omega\) to the discrete frequency \(\omega\).
• Crucially, the stable region (LHP in s-plane) should map to the stable region (inside the unit circle in z-plane).

Step 1: Digital Filter Specifications

• Define the desired characteristics of the digital filter:
  • Critical Frequencies:
    • Passband edge: \(\omega_p\)
    • Stopband edge: \(\omega_s\)
    • (Normalized digital angular frequencies, e.g., \(0 \le \omega \le \pi\))
  • Tolerances:
    • Max. passband ripple/attenuation: \(A_p\) (dB) or \(\delta_p\)
    • Min. stopband attenuation: \(A_s\) (dB) or \(\delta_s\)

Step 2: Method 1 - Impulse Invariance

• Concept: Match the impulse response: \(h_{digital}[n] \approx T \cdot h_{analog}(nT)\).
• Mapping: Maps s-plane pole \(s_k\) to z-plane pole \(z_k = e^{s_k T}\).
  • Uses partial fraction expansion of \(H(s)\).
• Pros: Preserves impulse response shape.
• Cons: Suffers from frequency aliasing if \(H(j\Omega)\) is not bandlimited below Nyquist frequency (\(F_s/2\)). Best for low-pass/band-pass filters.

Step 2: Method 2 - Bilinear Transformation

• Concept: An algebraic substitution mapping the \(j\Omega\) axis onto the unit circle.
• Mapping Formula: \[s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\] (T = sampling period, often set to 2 for simplicity during derivation).
• Pros:
  • No aliasing: Maps entire \(j\Omega\) axis to the unit circle (\(-\pi \le \omega \le \pi\)).
  • Preserves stability: Maps LHP (stable s-region) to inside the unit circle (stable z-region).
• Cons: Introduces non-linear frequency warping.

Understanding Frequency Warping (Bilinear)

• The relationship between analog (\(\Omega\)) and digital (\(\omega\)) frequencies is non-linear: \[\Omega = \frac{2}{T} \tan\left(\frac{\omega}{2}\right)\] or equivalently: \[\omega = 2 \arctan\left(\frac{\Omega T}{2}\right)\]
• This compresses the infinite analog frequency axis onto the finite digital frequency range \([-\pi, \pi]\).
• The mapping is non-uniform, most distorted near \(\omega = \pm \pi\).

Step 3: Frequency Pre-warping (Bilinear Only)

• Problem: Frequency warping distorts the locations of \(\omega_p\) and \(\omega_s\).
• Solution: Pre-warp the digital specifications (\(\omega_p, \omega_s\)) into corresponding analog frequencies (\(\Omega_p, \Omega_s\)) before designing the analog filter.
• Formula: Use the warping relationship: \[\Omega_p = \frac{2}{T} \tan\left(\frac{\omega_p}{2}\right)\] \[\Omega_s = \frac{2}{T} \tan\left(\frac{\omega_s}{2}\right)\]
• Use these \(\Omega_p, \Omega_s\) values for the analog design step.

Step 5: Apply Digital Transformation

• Convert the designed analog filter \(H_a(s)\) to the digital domain \(H(z)\):
  • If Impulse Invariance: Use partial fractions and the \(s_k \rightarrow e^{s_k T}\) mapping.
  • If Bilinear Transformation: Substitute \(s = \frac{2}{T} \frac{1 - z^{-1}}{1 + z^{-1}}\) into \(H_a(s)\) and simplify algebraically.
• Result: The digital filter transfer function: \[H(z) = \frac{\sum_{k=0}^{M} b_k z^{-k}}{1 + \sum_{k=1}^{N} a_k z^{-k}}\]

Step 6: Realization (Implementation)

• The transfer function \(H(z)\) corresponds directly to a difference equation: \[y[n] = \sum_{k=0}^{M} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]\]
• \(x[n]\): Input sequence
• \(y[n]\): Output sequence
• \(b_k, a_k\): Filter coefficients derived from \(H(z)\).
• The feedback term (sum involving \(y[n-k]\)) makes the impulse response potentially infinite (IIR).

Step 0: The Analog Domain - Laplace Transform

• Context: Continuous-time signals and Linear Time-Invariant (LTI) systems.
• Tool: The Laplace Transform converts differential equations (time domain) to algebraic equations (s-domain).
• Transfer Function \(H(s)\): Represents the system in the s-domain. \[H(s) = \frac{Y(s)}{X(s)}\] where \(s = \sigma + j\Omega\) is the complex frequency (\(\Omega\) = analog angular frequency).

Step 0: The Analog Domain - s-Plane Analysis

• s-Plane: Complex plane where \(H(s)\) is analyzed.
• Poles & Zeros: Roots of the denominator & numerator of \(H(s)\). Determine filter behavior.
• Frequency Response: \(H(j\Omega)\).
• Stability: All poles must be in the Left-Half Plane (LHP) (\(\text{Re}\{s\} < 0\)).

Step 0: The Analog Domain - Filter Prototypes

• Standardized analog filter approximations:
  • Butterworth: Maximally flat passband.
  • Chebyshev I/II: Equiripple passband/stopband.
  • Elliptic (Cauer): Equiripple passband & stopband (sharpest transition).
• These form the basis for many IIR digital filter designs.

The Bridge: Analog-to-Digital Mapping

• Core Task: Convert \(H(s)\) (s-plane) \(\rightarrow\) \(H(z)\) (z-plane).
• Relates continuous frequency \(\Omega\) to discrete frequency \(\omega\).
• Stable LHP (s-plane) \(\rightarrow\) Stable interior of unit circle (z-plane).
• Common Method: Bilinear Transformation (preserves stability, avoids aliasing, but warps frequency).

Step 1: Digital Filter Specifications

• Define the desired characteristics of the digital filter:
  • Sampling Frequency: \(F_s\) (Hz)
  • Critical Frequencies (Digital):
    • Passband edge: \(f_p\) (Hz) \(\implies \omega_p = 2\pi f_p / F_s\) (radians/sample)
    • Stopband edge: \(f_s\) (Hz) \(\implies \omega_s = 2\pi f_s / F_s\) (radians/sample)
    • Note: Scipy often uses frequencies normalized to Nyquist: \(W_p = f_p / (F_s/2)\), \(W_s = f_s / (F_s/2)\)
  • Tolerances:
    • Max. passband ripple/attenuation: \(A_p\) or gpass (dB)
    • Min. stopband attenuation: \(A_s\) or gstop (dB)

Step 1: Example Specifications (Python)

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# Filter Specifications
fs = 10000  # Sampling frequency (Hz)
fp = 1000   # Passband edge frequency (Hz)
fs_stop = 1500 # Stopband edge frequency (Hz)
gpass = 1    # Passband ripple (dB) - max loss
gstop = 40   # Stopband attenuation (dB) - min attenuation

# Normalize frequencies to Nyquist frequency (fs/2)
wp = fp / (fs / 2)
ws = fs_stop / (fs / 2)

print(f"Sampling Frequency: {fs} Hz")
print(f"Normalized Passband Edge: {wp:.3f} (pi rad/sample)")
print(f"Normalized Stopband Edge: {ws:.3f} (pi rad/sample)")
print(f"Max Passband Loss: {gpass} dB")
print(f"Min Stopband Attenuation: {gstop} dB")

Steps 2-4: Finding Order & Designing (Conceptual)

• Transformation: Bilinear Transform is implicitly used by many scipy.signal functions for IIR design.
• Frequency Pre-warping: Handled internally by Scipy when given digital frequency specs.
• Analog Design: Scipy calculates the required analog prototype order and parameters based on the (internally pre-warped) specs.
• Goal: Find filter order \(N\) and digital cutoff frequency \(W_n\).

Steps 2-4: Finding Order & Wn (Python - Butterworth Example)

• Use scipy.signal.buttord to find the minimum order and cutoff frequency for a digital Butterworth filter meeting the specs.

# --- Continued from previous code block ---

# Find the minimum order and cutoff frequency
N, Wn = signal.buttord(wp, ws, gpass, gstop, analog=False)

print(f"\nMinimum Filter Order (N): {N}")
# Wn is the digital cutoff freq (normalized) for butter function
print(f"Butterworth Natural Frequency (Wn): {Wn:.3f} (pi rad/sample)")

• analog=False specifies we want parameters for a digital filter.
• Wn is the digital cutoff frequency (scalar for LP/HP, tuple for BP/BS) required by the butter function. It’s often close to wp.

Step 5: Designing the Digital Filter (Python - Butterworth)

• Use scipy.signal.butter (or cheby1, cheby2, ellip) to get the filter coefficients \(b\) (numerator) and \(a\) (denominator) of \(H(z)\).

# --- Continued from previous code block ---

# Design the digital Butterworth filter
# btype='low' for lowpass
b, a = signal.butter(N, Wn, btype='low', analog=False, output='ba')

print("\nFilter Coefficients:")
print(f"Numerator (b): {np.round(b, 5)}")
print(f"Denominator (a): {np.round(a, 5)}")
# H(z) = (b[0] + b[1]z^-1 + ...)/(a[0] + a[1]z^-1 + ...)
# Note: a[0] is usually 1 after normalization

• output='ba' gives numerator/denominator coefficients. Other options: 'sos' (second-order sections - numerically better), 'zpk' (zeros, poles, gain).

Visualizing Frequency Response (Python)

• Calculate the frequency response using scipy.signal.freqz.
• Plot magnitude (in dB) and phase.

# --- Continued from previous code block ---

# Calculate frequency response
w, h = signal.freqz(b, a, worN=8000, fs=fs) # Use fs for physical freq axis

# Plot Magnitude Response
plt.figure(figsize=(10, 6))
plt.subplot(2, 1, 1)
plt.plot(w, 20 * np.log10(abs(h)))
plt.title('Digital Butterworth Filter Frequency Response')
plt.ylabel('Amplitude [dB]')
plt.grid(True)
plt.axvline(fp, color='green', linestyle='--', label=f'fp={fp}Hz') # Passband edge
plt.axvline(fs_stop, color='red', linestyle='--', label=f'fs={fs_stop}Hz') # Stopband edge
plt.axhline(-gpass, color='green', linestyle=':', label=f'-{gpass}dB')
plt.axhline(-gstop, color='red', linestyle=':', label=f'-{gstop}dB')
plt.legend()
plt.xlim(0, fs/2)
plt.ylim(-60, 5) # Adjust ylim as needed

# Plot Phase Response
plt.subplot(2, 1, 2)
angles = np.unwrap(np.angle(h))
plt.plot(w, np.degrees(angles))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Angle [degrees]')
plt.grid(True)
plt.xlim(0, fs / 2)

plt.tight_layout()
plt.show()

Step 6: Realization & Filtering (Python Example)

• The coefficients b, a define the difference equation.
• Use scipy.signal.lfilter to apply the filter to a signal.

# --- Create a sample signal ---
T = 1.0 # seconds
n_samples = int(T * fs)
t = np.linspace(0, T, n_samples, endpoint=False)

# Signal: Sine wave inside passband + sine wave inside stopband
sig_pass = 0.8 * np.sin(2 * np.pi * fp/2 * t) # 500 Hz
sig_stop = 0.4 * np.sin(2 * np.pi * (fs_stop + 500) * t) # 2000 Hz
signal_in = sig_pass + sig_stop

# --- Apply the filter ---
filtered_signal = signal.lfilter(b, a, signal_in)

# --- Plot signals ---
plt.figure(figsize=(10, 5))
plt.plot(t, signal_in, label='Input Signal', alpha=0.7)
plt.plot(t, filtered_signal, label='Filtered Signal', linewidth=2)
plt.title('Filtering Example')
plt.xlabel('Time [sec]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)
plt.ylim(-1.5, 1.5)
plt.show()

# Observe: The 500Hz component passes, the 2000Hz component is attenuated.